function [varargout] = lastdimfun(func, varargin)
%LASTDIMFUN Apply function to each subarray along the last dimension of array

% make sure for all input arrays:
% 1. number of dimensions greater than 1
% 2. sizes of the last dimension are equal
for iArgIn=1:nargin-1
    sz = size(varargin{iArgIn});
    assert(length(sz) > 1);
    if iArgIn == 1
        nLastDim = sz(end);
    else
        assert(sz(end) == nLastDim);
    end
end

% select the i-th subarray of each input array and pass them to func
for iLastDim = 1:nLastDim
    funcArgIn = cellfun(@(A)selectsubarr(A, iLastDim), varargin, 'UniformOutput', false);
    if nargout > 0
        [funcArgOut{1:nargout, iLastDim}] = func(funcArgIn{:});
    end
end

% convert each output of func to an array
for iArgOut = 1:nargout
    varargout{iArgOut} = cell2mat(funcArgOut(iArgOut, :));
end
end

function [ASub] = selectsubarr(A, i)
ind = totalindofdims(A); % TODO: optimize for speed
ind{end} = i;
ASub = A(ind{:});
end